New proposal for numerical simulations of ^-vacuum like systems 
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We propose a new approach to perform numerical simulations of 0-vacuum like systems, test it 
in two analytically solvable models, and apply it to CP 3 . The main new ingredient in our approach 
is the method used to compute the probability distribution function of the topological charge at 
6 = 0. We do not get unphysical phase transitions (flattening behavior of the free energy density) 
and reproduce the exact analytical results for the order parameter in the whole 6>-range within a 
few percent. 
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Quantum Field Theories with a topological term in 
the action are a subject of interest in high energy parti- 
cle physics and in solid state physics since a long time. 
In particle physics, these models describe particle inter- 
actions with a CP violating term. The extremely small 
experimental bound for the CP violating effects in QCD 
(strong CP problem) is still waiting for a convincing the- 
oretical explanation 0. In solid state physics, chains of 
half-integer quantum spins with antiferromagnetic inter- 
actions arc related to the two-dimensional 0(3) nonlinear 
sigma model with topological term at 9 ~ tt. It has been 
argued that this model presents a second order phase 
transition at 9 = tt, keeping its ground state CP sym- 
metric (Haldane conjecture) 

Non-perturbative studies of field theories with a 9- 
vacuum term are enormously delayed because of the com- 
plex character of the euclidean action which forbids the 
application of all standard Monte-Carlo algorithms. Be- 
sides this, lattice QCD lacks from a simple consistent 
definition of topological charge. 

The partition function Zy{9) of any 9- vacuum like 
model in a finite space-time lattice volume V can be 
written, up to a normalization constant, as the discrete 
Fourier transform of the probability distribution function 
(p.d.f.) of the topological charge at 9 = 0: 

Z V (9) = Y,Pv(n)e Wn , (1) 

n 

where pv (n) is the probability of the topological sector 
n and the sum runs over all integers n. In almost all 
practical cases the sum in (|l|) has a number of terms 
of order V since the maximum value of the topological 
charge at finite volume is of this order. 

Since efficient algorithms for numerical simulations of 
physical systems with complex actions are not yet avail- 
able, the only a priori reliable numerical scheme to an- 
alyze the thermodynamics of ^-vacuum like models goes 
through the determination of the p.d.f. of the topological 
charge, pv{n), and the evaluation of its Fourier trans- 



form (0). But this is a difficult task due to the following 
two technical reasons, which will be clarified later: i) 
any numerical determination of pv (n) suffers from sta- 
tistical fluctuations j|, and small errors in pv(n) can 
induce enormous relative errors in the determination of 
a quantity as Zy (9) which is an extremely small number 
of order e~ v , ii) even if we were able to evaluate py(n) 
with infinite precision, the sum in (|l|) contains terms that 
differ by many orders of magnitude, running from 1 to 
e~ v §. 

In few specific cases one can overcome the sign prob- 
lem ||. However previous attempts by other groups to 
simulate 9- vacuum like systems ^, |(| , were based on the 
numerical determination of the p.d.f. of the topological 
charge straightforwardly, by standard simulations, or by 
more sophisticated methods based on the use of multi- 
binning and re-weighting techniques. In all these at- 
tempts, artificial phase transitions at a 9 C decreasing with 
the lattice volume were observed for the two-dimensional 
U(l) gauge theory at strong coupling as well as CP^ 
models. The origin of this artificial behavior, which fol- 
lows from a flattening behavior of the free energy for 
9- values larger than a certain 9 C , was analyzed in 0. 
Both groups agreed that the observed behavior was pro- 
duced by the small statistical errors in the determination 
of the p.d.f. of the topological charge, the effect of which 
became more and more relevant as the lattice volume 
was increased. In Ref. M it was also noticed that by 
smoothing the p.d.f. flattening disappears. 

The purpose of this paper is to introduce a new nu- 
merical approach to simulate 9- vacuum like models. This 
approach is based on a new method to compute the p.d.f. 
of the topological charge and the use of a multi-precision 
algorithm in order to compute the sum in ([!]) with a pre- 
cision as high as desired. 

For reasons which will become apparent in what fol- 
lows, let us write the partition function (|l|) as a sum 
over the density of topological charge x n = n/V and set 
Pv{ n ) — cx p[~ V fv{xn)], where is a smooth inter- 
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polation of — \/V lnpv(x n ): 

Z v {9) = J2 e' y f v ^e WVx -. (2) 

Equation (||) defines a 2ir periodic function of 9. Since 
CP is a symmetry of the action at 9 — and 9 = ir, 
fv(x) will be an even function. We will assume that CP 
is realized in the vacuum at 9 = since otherwise the 
theory would be ill-defined at 9 ^ ||. This implies 
that exp[— Vfy(x n )] will approach a delta distribution 
centered at the origin in the infinite volume limit. In 
some exceptional cases, as QCD in the chiral limit, the 
function /y (x) is not defined since any topological sector 
with non vanishing charge has a vanishing probability. 
However, this is a trivial case in which the theory is in- 
dependent of 9. 

Let us consider the partition function (||) in the com- 
plex 9 plane, in particular on the imaginary axis 9 = —ih, 
and let f{x) be the infinite volume limit of /y(x n ). All 
the terms entering equation ^ are positive definite, and 
then in the infinite volume limit the free energy is given 
by the saddle point. Assuming that f(x) has first deriva- 
tive for any x except at most in isolated points, we can 
write the saddle point equation: 

f{x) = h (3) 

which gives the external " magnetic field" h as a function 
of the density of topological charge x. 

Our proposal to compute the function f(x) is based in 
the following three steps: 

i. To perform standard numerical simulations of our 
system at imaginary 9 = —ih and to measure the 
mean value of the density of topological charge as 
a function of h with high accuracy (tipically a frac- 
tion of percent). This is feasible since the system 
we have to simulate has a real action. Then, Eq. (||) 
is used to get a numerical evaluation of f'{x). 

ii. To get f(x) we have to integrate f'(x). Between 
the possible ways to do this integral, we decided to 
fit f'(x) by the ratio of two polynomials, whose or- 
der is chosen to obtain a high quality fit, and then 
to perform analytically the integral of the fitting 
function. In this way we get a very precise deter- 
mination of f{x), which allows us to compute the 
p.d.f. in a range varying several thousands orders of 
magnitude. This is the main advantage of our ap- 
proach when compared with other methods based 
on a direct computation of pv(n). 

iii. To use a multi-precision algorithm to compute the 
partition function (^) using as input the function 
fix) previously determined. 

The function f(x) obtained in step ii) suffers from sta- 
tistical and systematic errors, the last coming from the 




FIG. 1: Magnetization versus the imaginary magnetic field 
in the one- dimensional antiferromagnetic Ising model at F = 
— 1/2 on a chain of 1000 sites: exact (dashed) and numerical 
(continuous) results. Statistical errors are smaller than 2%. 

fact that the saddle point equation (0) has finite volume 
corrections. An analysis of these errors for the models 
and sizes we have studied (see below), shows that sys- 
tematic errors due to finite volume effects are smaller 
than the statistical ones in the whole relevant range of x. 
This is the reason why we will replace fv(x n ) in equation 
(Q) by its asymptotic value f(x n ) in what follows. This 
substitution has no effect in the infinite volume limit at 
imaginary 9 and we are assuming that the same holds for 
real 9. 

Before presenting our results for the various testing 
models we have analyzed, let us briefly discuss how errors 
in the determination of f(x) can propagate to Zy(9). 
This is an important point since, as pointed out in 
the artificial phase transitions found in U(l) and CP 
were caused by the statistical errors in the determination 
of the p.d.f. 

To this end, let /y(x n ) be the exact value at a given 
lattice volume V of the function which parameterizes the 
p.d.f. and be A/y(x n ) a given deviation from the ex- 
act value. If we denote by Z'y(9) the partition function 
computed with fv{x n ) + A/y(x„), a simple calculation 
tells us that this partition function is related with the 
exact partition function Zy{9) as follows: 

Z' v {9) = Z v (9)(e- VAfv ^ ) ) . (4) 

We can at this point analyze two extreme cases. First 
let us assume that Afv(x n ) vanishes everywhere ex- 
cept at a given x m . Taking into account that the par- 
tition function Zy{9) which enters in the denomina- 
tor of the expectation value {er v ' v ^ n ') should be- 
have as e~ Vgv ( e \ where gv(@) is the free energy density 
(#v(0) = 0), we should get 

{e -VA fv ( Xn)) = l 

+ 2e - y (/v-(v)-sv(e))( e -VA.M^) _ l) C os(m6»). (5) 
Since gv(0) is an increasing function of 9, it will be 
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FIG. 2: Topological charge versus 6 in the two-dimensional 
U{1) model at = and (3 = 0.6 on a 80 x 80 lattice. Sta- 
tistical errors are not visible at this scale. The exact result 
(dashed curve) cannot be distinguished from the numerical 
result (continuous curve). 



smaller than fv(x m ) near the origin, and therefore the 
free energy density computed with the modified partition 
function will differ from the exact one of a small quantity 
of order y in this region. However at larger values of 9 
the function gv{9) can become larger than fv{x m ) and 
in such a case the correction to the free energy density 
will be finite or, even worst, the modified partition func- 
tion can become negative. The other extreme case is that 
in which we assume that the error Afv(x n ) is constant. 
Under such assumption Eq. implies that the error in 
the free energy density gv{9) will also be Afv- 

The previous discussion, as the results of sug- 
gests that correlated errors propagate in a less dramatic 
way than uncorrelated ones, and this would be a good 
scenario for numerical methods which, as in our ap- 
proach, produce correlated errors in the determination 
of fv(xn)- We have checked this interesting issue in the 
two-dimensional U(l) model at strong coupling (see be- 
low) and verified that in fact correlated errors in fv{x n ) 
induce errors of the same order in the free energy density 

9v(0)- 

To test these ideas we have analyzed three models: the 
one-dimensional antiferromagnetic Ising model within an 
external imaginary magnetic field, the two-dimensional 
compact U(l) model with topological charge, and CP 3 in 
two dimensions. The coupling to the imaginary magnetic 
field in the Ising model can be written as i9ksT/2 ^ Si, 
where ks is the Boltzmann constant and T the tem- 
perature. For an even number of spins, the quantity 
1/2 Si is an integer taking all values between —N/2 
and N/2, and therefore it can be seen as a quantized 
charge. Furthermore, the theory has a Zi symmetry at 
9 = and 9 = it which is the analogous of CP in field the- 
ory We use the notation F = J/ksT, where J is the cou- 
pling constant between nearest neighbors. The transfer 
matrix technique allows to solve analytically the model. 



FIG. 3: Topological charge versus 6 in the two-dimensional 
U(l) model at (5 = 0. The continuous curve is the exact 
result, the dashed curve is the result obtained by substituting 
f(x) by f(x)(l + i sin(a; 2 )), and the dotted curve is the result 
obtained by adding a random error of the order of 0.1% to 
/(*)• 

For antiferromagnetic couplings, F < 0, the magnetiza- 
tion is an analytic function of 9 between —tt and tt. At 
9 = tt the system shows a first order phase transition with 
a non-vanishing magnetization. From a numerical point 
of view the determination of the free energy density and 
order parameter through equation (|]J) in this model has 
the same level of complexity of more sophisticated mod- 
els. Furthermore, in contrast to two-dimensional U(l) 
gauge theory, where the p.d.f. of the topological charge 
is nearly gaussian, the non-gaussian behavior of the p.d.f. 
of the mean magnetization in the antiferromagnetic Ising 
model makes this model a good laboratory to check the 
reliability of our approach. Fig. [I] shows our numerical 
results for the order parameter versus 9 for a linear chain 
of 1000 spins and F = —1/2. Statistical errors were esti- 
mated by doing 10 samples of the numerical results and 
applying a Jack-Knife analysis. The p.d.f. of the order 
parameter for such a system takes values in a range of 
around 2000 orders of magnitude. Notwithstanding that, 
we are able to reproduce the order parameter in the whole 
9 interval within a few percent. 

The two-dimensional compact U(l) gauge model with 
0-angle at strong coupling constitutes another interesting 
check because we can compare the goodness of our ap- 
proach with the other existing simulations which showed 
artificial behavior with a fictitious phase transition mov- 
ing to the origin when increasing the lattice volume. 
Fig. ^ displays our results for the topological charge den- 
sity versus 9 in a 80x80 lattice at j3 = and (3 = 0.6. 
We are able to reproduce the exact result within a few 
per thousand in the whole 9 interval. The agreement 
between analytical and numerical results is actually im- 
pressive. Furthermore the flattening found in || for the 
free energy density in relatively small lattices is absent 
in our simulations even in the 80x80 lattice. 

To test how different kind of errors in the determina- 
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FIG. 4: Topological charge versus 6 in the two-dimensional 
CP 3 model at = 0.6 on a 100 x 100 lattice. Continuous line 
(discontinuous line) reports the results obtained fitting f'(x) 
with a polynomial (the ratio of two polynomials). Statistical 
errors are at 1% level. 



tion of the function f(x) which defines the p.d.f. of the 
density of topological charge can affect the determination 
of the free energy and order parameter, we have added to 
the measured f{x) a random relative error of order 10~ 3 . 
Fig. U shows the order parameter obtained in this way. 
As can be seen a small but random error in f(x) propa- 
gates to the order parameter in a very dramatic way and 
makes the calculation meaningless. Contrary to that if, 
in order to simulate a correlated relative error of order 
up to 50%, we replace the measured f(x) by the (even) 
function f(x) (l + 0.5 sin(x 2 )) , the result for the order 
parameter is practically indistinguishable from the exact 
value for 9 < 7r/2, and the maximum deviation is about 
25%, at 9 = 7r (see Fig. ||). We conclude that random 
errors in f(x) propagate in a very dramatic way but cor- 
related errors do not, and this helps one to understand 
why our approach works so well. 

The last model we have analyzed is CP 3 in two dimen- 
sional euclidean space. It is the standard wisdom that 
this model shares many qualitative features with QCD. 
Even if it has not been analytically solved we believe it is 
worthwhile to compare our results with previous existing 
numerical simulations. We studied the lattice formula- 
tion that makes use of an auxiliary U(l) field. Also in 
this model, the previous numerical simulations gave ar- 
tificial phase transitions with a flattening behavior for 
the free energy density at a 9 C decreasing with the lattice 
volumei, §. 

Fig. [| shows our results for the order parameter versus 
9 at /3 = 0.6 on a 100 2 lattice. We have chosen this par- 
ticular [3 value in order to compare directly with results 
reported in || . In our simulations we have no trace of the 
fictitious phase transition found in || . Furthermore the 
order parameter is clearly different from zero at 9 = 7r, 
hen ce CP is spontaneously broken at this value. An 



eluding a study of the behavior of the order parameter 
in the continuum limit, and an analysis of statistical and 
systematic errors involved in our approach, will be pub- 
lished elsewhere. What is interesting to notice here is 
that our results do not suffer from artificial phase tran- 
sitions, caused by statistical errors in the determination 
of the p.d.f. 

In the three models analyzed the finite size effects in 
f'(x) cannot be appreciated since they are completely 
masked by the small statistical errors. For instance, finite 
size effects can be exactly computed in the Ising model: 
they are exponentially small with the lattice size. This is 
a general feature of non-critical systems. However, vol- 
ume effects might be troublesome in the analysis of the 
continuum limit. Concerning systematic errors due to the 
choice of a particular fitting function for f'(x), the differ- 
ence between the numerical and exact results for the Ising 
and compact U(l) models (beside the statistical errors) 
reported in Figs. 1 and 2 give us an idea on the order 
of magnitude of these errors. Of course systematic errors 
can depend on the model as well as on the parameters. 
In CP 3 at = 0.6 we did also a 5-parameters polyno- 
mial fit of the data for f'{x) in the relevant x-interval. 
The discrepancy between the topological charge density 
obtained in the two cases is at most 7% (see Fig. 4). 

Similar ideas to those presented in this work have been 
proposed and promisingly applied to a matrix model of 
QCD at finite density in |lC|| . 
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opc ji question is how CP is realized in the continuum 
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